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I.  SUMMARY 


The  research  performed  under  this  contract,  during  the  period  of  12  January  1984 
through  11  January  1985,  involved  several  tasks  including; 

(a)  Development  of  computer  programs  for  an  elastic  Kircholf  method  to  handle  SH,  SV, 
and  P  interactions. 

(b)  Simulation  of  complex  structures  for  given  velocity  models  and  source  properties  and 
testing  against  numerical  methods. 

(c)  Development  of  analytical  theory  to  treat  high-frequency  propagation  in  laterally 
varying  structure. 

This  report  will  discuss  two  aspects  of  our  progress,  namely  the  analytical  develop¬ 
ments  in  treating  laterally  varying  structure  and  numerical  extensions  in  theory  to  han¬ 
dle  the  air-ground  interface  problem. 

In  section  II,  we  present  some  extensions  of  WKBJ  and  generalized  ray  theory  for 
line  and  point  sources,  Helmberger  et.  al  (1985).  In  particular,  some  changes  in  existing 
methods  are  required  to  construct  synthetics  for  broad-band  signals  in  many  situations, 
especially  when  energy  can  reach  the  receiver  by  up-going  as  well  as  down-going  paths. 
This  can  be  accomplished  by  allowing  locally  dipping  structure  and  making  some 
rnodilications  to  generalized  ray  theory.  Local  ray  parameters  are  expressed  in  terms  of 
a  global  reference  which  allows  a  de  Hoop  contour  to  be  constructed  for  each  generalized 
ray  with  the  usual  application  of  the  Cagniard-de  Hoop  technique.  Several  useful 
approximations  of  ray  expansions  and  WKBJ  theory  are  presented.  Comparisons  of  the 
synthetics  produced  by  these  two  basic  techniques  alone  or  in  combination  with  known 
solutions  demonstrates  their  reliability  and  limitations.  In  section  III,  we  present  a  tech¬ 
nique  for  improving  the  stability  of  finite-  dilTcrence  calculations  at  a  free  surface,  Vidale 


.  *>  . 


and  Clayton  (1985)  where  a  stable  free-surface  boundary  condition  is  derived  for  finite- 
difference  solutions  of  the  two-dimensional  elastic  wave  equation.  Lateral  velocity  and 
density  variations  at  the  free  surface  are  correctly  treated.  The  method  is  implicit,  but 
only  requires  a  simple  pentadiagonal  system  solver  to  implement.  When  coupled  with 
second-  and  fourth-order  interior  solutions,  the  overall  problem  appears  to  be  stable  for 
shear  to  compressional  wave  velocity  greater  than  0.01  and  0.02,  respectively.  The 
method  is  compared  with  other  published  algorithms  for  accuracy  and  stability. 
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II.  Notes  On  Wave  Propagation  In  Laterally  Varying  Structure 

Considerable  progress  has  been  made  recently  in  speeding-up  the  synthesizing  of 
seismograms  with  the  introduction  of  WKBJ  and  Gaussian  beam  methods,  see  Chapman 
(1978)  and  Cerveny  et  al.  (1982).  These  methods  have  proven  highly  useful  in  generali¬ 
zations  to  laterally  varying  structure,  especially  at  high  frequency,  see  for  example 
Frazer  and  Phinney  (1980).  However,  in  the  construction  of  longer  periods  (long  period 
WW’SSN  seismograms)  we  many  times  are  interested  in  more  complete  solutions,  since 
the  beginning  portion  of  surface  waves  become  important,  see  Grand  and  Helmberger 
(1981a).  A  complete  set  of  ray  parameter  contributions  is  required  to  construct  seismo¬ 
grams  in  this  situation.  In  particular,  one  needs  to  consider  ray  paths  leaving  the  source 
horizontally,  a  case  where  the  WKBJ  method  breaks  down.  We  can  avoid  this  problem 
by  applying  a  mixture  of  generalized  ray  theory,  GRT,  and  WKBJ  or  Disk  rays  as 
defined  by  Wiggins  (1976). 

A  simple  example  of  this  procedure  is  given  in  Fig.  1  where  we  show  schematically 
how  to  construct  the  step  response  for  a  smooth  velocity  model  approximated  by  a  stack 
of  liomogeneous  layers.  We  suppose  that  a  velocity  model  can  be  chosen  such  that  the 
step  response  remains  a  step  at  all  receiver  positions.  The  simulation  of  this  step  can  be 
achieved  by  summing  the  response  from  three  energy  paths;  namely,  the  direct,  the 
reflected  from  just  below  the  source  or  reference  plane,  and  the  diving  WKBJ  contribu¬ 
tion.  All  three  paths  contain  a  product  of  the  transmission  coefficients  above  the  source. 
The  WKBJ  path  includes  the  transmission  coefficients  across  the  reference  plane,  taken 
as  the  interface  below  the  source.  We  have  included  a  diagram  of  the  0(/ )  vs.  p  curve 
in  Fig.  1  for  reference  as  it  clearly  shows  that  the  diving  path  contributes  little  except  at 
the  larger  distances.  At  the  nearest  distance,  position  1,  the  direct  ray  dominates.  The 
reflected  path  contributes  some  as  critical  angle  is  approached.  At  still  larger  distances, 
position  3,  a  head  wave  along  the  bottom  of  the  reference  interface  develops  followed  by 
the  critically  reflected  pulse.  The  head  wave  contribution  is  included  in  the  reflected 


Schematic  picture  of  ray  paths  and  summation  response  for  a  smootli  velo- 
;radient  specified  by  a  layered  stack.  Reflected  and  headwave  paths  for  case 
omitted. 
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response  since  it  is  associated  with  the  reflected  generalized  ray.  At  large  ranges,  the 
WKBJ  contribution  becomes  increasingly  dominant.  Note  that  the  WKBJ  response 
turns  off  at  the  same  time  as  the  head  wave  starts  because  the  transmission  coeiTicient 
drops  to  zero.  Essentially,  combining  the  generalized  rays  and  WKBJ  response  elim¬ 
inates  the  truncation  phase  and  avoids  the  turning  point  break-down  of  the  WKBJ 
theory. 

In  testing  the  accuracy  of  the  above  procedure  it  is  quite  useful  to  generate  the  step 
response  for  models  for  which  the  answer  is  known.  Thus,  we  begin  with  a  homogeneous 
fluid  whole  space  with  a  point  source  excitation  yielding  a  step  response  at  all  positions 
with  l/(distance)  decay.  We  next  impose  a  spherical  coordinate  system  with  many  thin 
shells  of  constant  velocity.  Applying  the  classical  earth-flattening  approximation  w'e 
obtain  a  model  with  a  smooth  velocity  increase  in  depth,  see  Helmberger  (1973).  The 
synthetics  generated  in  Fig.  2  are  from  such  a  model  with  the  exact  step  responses  indi¬ 
cated  by  the  dotted  lines  in  the  bottom  panel.  This  panel  also  displays  the  response 
after  summing  the  complete  set  of  generalized  rays;  direct  plus  rays  reflected  upward 
from  all  the  interfaces  below  the  source.  The  GRT  response  at  the  largest  distance  shows 
the  most  roughness  for  times  near  the  direct  arrival  when  the  interaction  with  the 
reflection  from  just  below  the  source  is  the  most  severe.  Similar  complexity  occurs  with 
the  hybrid  method  except  that  the  diving  energy  is  smoother  with  WKBJ.  Short  period 
synthetics  generated  from  these  step  responses  become  quite  dirty  and  simple  geometric 
ray  theory  yields  cleaner  results.  However,  for  most  studies  the  advantage  of  being  able 
to  include  the  radiation  pattern  appropriate  for  earthquake  sources,  or  shear  disloca¬ 
tions,  far  outweighs  the  disadvantage  of  the  noise  generated  by  the  hybrid  method.  For 
example,  consider  the  SH-radiation  from  a  dip-slip  event  where  the  up-going  radiation 
has  opposite  polarity  from  the  down-going  energy,  see  Helmberger  (1973).  In  short,  the 
sum  as  displayed  in  Fig.  2  becomes  more  interesting  when  the  direct  ray  trace  has  oppo¬ 
site  sign  from  the  other  two. 
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We  could  probably  improve  the  response  at  the  time  the  three  energy  paths  inter¬ 
fere  the  most  vigorously  by  including  a  few  more  GR’s  and/or  by  lowering  the  reference 
boundary  for  the  WKBJ  contribution.  However,  we  are  particularly  interested  in  more 
realistic  earth  models  with  a  sedimentary  cover  over  bedrock  or  a  crust-over-mantle 
structure  providing  natural  reference  boundaries.  Thus,  we  propose  using  GR’s  to  com¬ 
pute  the  start  of  the  Love  waves,  and  WKBJ  to  generate  the  responses  returning  from 
deeper  structure.  This  approach  proved  effective  in  studying  the  structure  and  evolution 
of  the  lithosphere  for  an  old  oceanic  plate,  Grand  and  Helmberger  (1984b).  It  would  be 
advantageous  to  treat  the  obvious  lateral  variation  encountered  in  such  studies. 
Although  the  real  world  is  truly  three  dimensional,  some  useful  progress  can  be  made  by 
examining  profiles  of  data  along  paths  of  symmetry  where  two-dimensional  idealizations 
are  appropriate.  We  will  address  such  models  in  this  paper. 

Our  strategy  is  similar  to  Wiggins  (1976)  in  that  we  will  use  a  combination  of  GRT 
and  DRT  to  generate  synthetics  and  justify  the  latter  by  demonstrated  accuracy. 

Review  of  Ray  interactions  with  nonplanar  structure 

Boundary  value  problems  involving  complicated  geometry  have  a  long  rather  unre¬ 
warding  history;  thus,  we  will  jump  directly  to  approximate  solutions  and  test  their  vali¬ 
dity  against  finite-difference  calculations  and  other  more  well-known  results.  Before 
addressing  the  dipping  layer  problem  it  is  instructive  to  examine  the  flat-layered  case 
and  emphasize  the  geometric  interpretation  of  generalized  ray  theory.  This  proves  par¬ 
ticularly  useful  for  constructing  generalizations  to  more  complicated  situations  since  the 
most  progress  in  understanding  these  problems  is  at  high  frequency.  Both  line  and  point 
sources  will  be  discussed  since  the  former  is  easier  to  understand  theoretically  and  for 
testing  against  numerical  results,  while  the  latter  is  necessary  for  studying  the  Earth. 


The  solution  of  the  scalar  wave  equation  asauming  line  source  excitation  for  gen¬ 
eralized  rays  as  given  by  Gilbert  and  KnopofT  (1961)  is 

=  (1) 

where  tg  =  R  /a,R^  =  +  r~,  and  a  =  velocity. 

is  defined  as  the  displacement  potential  with  the  index  L  used  to  remind  our¬ 
selves  of  the  line  excitation.  A  high  frequency  approximation  of  (1)  is 

)/(<-<  (2) 

and  the  motion  decays  with  distance  as  the  \/W .  The  solution  to  the  interface  problem 
setup  displayed  in  Fig.  3a  is 

♦i(r,.,0  =  Im(-L-^r(p))  (3) 

tjj  at 

where 


t  =  p(di  -t-  cfj)  +  (^) 

=  (4  - 

T  (p)  =  Transmission  coefficient 

The  symbol  (Im)  indicates  the  imaginary  part  of  the  complex  product  of  the  func¬ 
tions  of  ray  parameter,  see  Helmberger  (1983)  for  example.  The  ray  parameter  appropri¬ 
ate  for  the  direct  arrival  path,  p,  ,  can  be  obtained  by 


-T'^Po)  =  0>“d  di 
dp 


=  d, - 


But  with 


sin©i  sin02 


(6) 


and,  therefore, 


Figure  3:  Diagram  displaying  the  geometric  spreading  of  ray  tubes  in  two  and  three 
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Ql  02 


we  see  that  the  ray  goes  from  the  source  to  the  receiver.  And 


to  —  R  i/oi  +  R  2/  0-2 

For  times  greater  than  ,  we  must  solve  t  for  complex  p  such  that  the  imaginary  parts 
of  pdi  and  etc.  cancel. 

The  behavior  near  p,  can  be  approximated  by  noting  that 


and  solving  for 


Thus, 


Note  that  from  (4) 


ip  -p.f -SCI  -<,)/(44). 


d^t  _  1 _ ^2 

rfp*  ViQi 


It  is  convenient  to  condense  the  various  factors  containing  p,  into 


Sl(p.)^—^ - 

I  y4  I 

dp 

which  we  call  the  spreading  factor.  Thus, 
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We  note  that  by  differentiating  Snell’s  law  we  obtain 


- L^e,  - - irfOj  . 

0|  an 


Substituting  this  expression  into  Si  'Mt  obtain 


S,  = 


r  "  ■  ■  ■  '  ■  ,  -  11  ■  1 1  . 

/?  0j  +  R  jd  ©2 


If  /?2  =  0>  we  obtain  the  whole  space  spreading  again  where  {RidQi)  is  just  the  width 
of  the  ray  tube  described  in  Fig.  3a.  A  correction  for  the  change  in  direction  is  required 
as  the  tube  crosses  the  interface,  namely 

(cos0|/cos©2)  . 

Thus,  the  denominator  of  (12)  is  again  the  width  of  the  ray  tube  at  the  receiver,  m 
Fig.  3a.  Substituting  into  (3)  we  obtain 

<^1  a  5,//(<-tJ/?c{r(p.))/(f-lJ>/2  (13) 

where  R^  indicates  the  real  part  operator. 


The  point  source  solution  for  the  same  problem  setup.  Fig.  3b,  is 

and  applying  the  same  first  motion  approximation  we  obtain  a  slightly  more  complicated 


spreading  factor  namely. 


V  r  7,  dp- 


/?jsin0i  R^itiQn 
sin©! 


and  note  that  by  letting  R  2=^  0  vie  obtain 


Sp  =  \/R 


(17) 


In  terms  of  area,  we  note  that 


Sp  = 


sinO]d0](f  4 


(i?isin@id^  +  /?2sin62d4)(/?id0|  4-  /f2</02 


COS01 

COS0O 


1/2 


(18) 


which  can  be  interpreted  as  the  incremental  element  of  area  at  the  source  divided  by  the 
projected  area  at  the  receiver  or  simply, 


Thus,  the  first  motion  behavior  becomes 


(19) 


<l»,(ritiO  =  ^ 


TT  -TtT 


5;.Re(r(pJ) 


(20) 


=  5,/f(/-<,)Re(r(pJ) 


More  complicated  solutions  to  multi-layered  models  in  terms  of  ray  summations  will  be 
discussed  later. 


Locally  Dipping  Structure 

Although  GRT  for  parallel  interfaces  has  been  well  developed  the  modifications  for 
nonplanar  structure  or  smoothly  varying  interfaces  has  not.  Some  of  the  difficulties 
encountered  for  the  simple  wedge  problem  have  been  discussed  by  Hudson  (1963).  Hong 
and  Helmberger  (1977)  constructed  a  solution  in  terms  of  generalized  rays  for  this  prob¬ 
lem  and  defined  a  method  of  ray  path  construction  compatible  with  the  usual  Cagniard- 
de  Hoop  formalism.  We  will  consider  the  direct  arrival  interacting  with  two  dipping 
interfaces  as  an  example  application.  The  problem  setup  is  displayed  in  Figure  -la  with 
the  response  given  by 


I 


I  (Pi)  rasVPa)!  dp 

”  1  - Z - 


where  pj  and  P2  are  defined  by  the  local  ray  parameter,  namely 


sin6i  sin62 

Pi  - - ,  P2— - 

Ql  02 


and  are  no  longer  equal.  However, 


where 


00  =  02  “t" 


with  0,  defining  the  change  of  the  slope  of  interface  (1)  relative  to  the  previous  refer 
ence  at  0i.  Performing  the  derivatives  discussed  in  the  previous  section  we  obtain 
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which  is  similar  to  (12)  and  has  the  same  interpretation.  The  travel  time  is  defined  by 


t  =  'EiPidi  Vi  hi ) 


with  the  definitions  of  rf,  and  />,-  given  in  Fig.  -la  as  the  projection  of  the  geonietrir  path 
onto  the  local  Cartesian  coordinates.  The  arrival  lime  can  be  determined  a.s  before  with 


defining 


Pi  =  P«,.P2  =  Po, 


etc.  Thus. 
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with 


Pm 


»‘“®m  ^  _  COS0„ 

— -  ana  9m  - - 


dt 


and  the  — - —  0  condition  leads  to  a  ray  going  from  the  source  to  the  receiver.  The 

dPm 


first-motion  approximation  becomes 


(24) 


Spreading  for  the  point  source  solution  becomes  slightly  more  complicated  than  in 
the  flat  ca.se,  but  allowing 

t/2 


(25) 


results  in  Sp  defined  by 


5p 


=  /lJl.  I  Jl 

V  r  I,,  '  rfp 


I  -1/2 


(26) 


reducing  to 


{A, /A  )‘/2  , 


at  the  direct  arrival  lime.  The  details  of  this  result  have  been  given  previously  by  Hong 
and  Melmberger  (197S).  Thus,  the  point  source  solution  for  the  geometry  given  in  Fig. 
4b  becomes 


♦'-7 


1 

TT 


r*""  I 


(27) 


where 


n(p)—  Ti2(p|)r23(p2) 


Numerical  evaluation  of  (27)  yields  the  geometric  result  but,  also,  retains  longer  period 
information  since  (dp/dt)  can  be  evaluated  along  the  de  Hoop  contour  in  the  usual 
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manner.  The  accuracy  of  constructing  broad-band  synthetics  applying  this  procedure  is 
well  known,  see  for  example  Apsel  and  Luco  (1983)  or  Burdick  and  Orcutt  (1978). 


and  the  summation  of  n  rays  is  required.  The  various  symbols  are  defined  below: 

»;(r  ,z  ,Q,t )  =  displacement  on  free  surface 

A/j  =  moment 

Pg  =  density 

D(t)  =  dislocation  history 

D  (t)  =  far-field  time  function 

A  i(©,X,6)  =  co620cosXsinj  -  l/2sin20sinX«tn  2^, 

A  2(0, X, 5)  =  -sin  ©cosXcosA  -  cos©sinX cos  28, 

©  =  strike  from  the  end  of  the  fault  plane 
X  =  rake  angle 
A  —  dip  angle 

r  =  distance  between  source  and  receiver 


16 


-  13- 


p  =  ray  parameter 


A  =  epicentral  distance  in  radians 

- I  =  correction  for  earth  flattening 

sinA  I 


6  =  shear  velocity 

and  where  the  vertical  radiation  patterns  are  given  by 


SH 


SIIo 


J_  IL 
0^  P 


+1  >  A 

*  -1  z  <  h 


n  = 


1/2 


The  correction  for  point  source  spreading  is  defined  by 

,-»/2 


(31) 


This  solution  is  similar  to  the  flat  case  and  we  can,  obviously,  construct  the  diving  ray 
response  for  a  smoothly  varying  structure  by  summing  the  primary  rays  as  discussed  in 
Fig.  2.  We  can  then  use  this  result  to  check  the  disk  ray  solution  which  can  be  obtained 
by  replacing  (30)  by 

=  SHi{p)^{p)Y.{^p  /6l)  (32) 

where  the  sum  is  taken  over  the  p  (t)  curve  as  described  by  Wiggins  (1976). 

For  a  simple  turning  ray  problem 


I  i£|  _ \ _ 

8t  |'•-r(p)| 


(33) 


where  r(p)  is  distance  reached  by  a  ray  defined  by  p,  see  Fig.  5.  Substituting  (33)  into 
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receiver 


Figure  5;  ftiiy  geometry  where  the  source  and  receiver  are  separated  by  r  and  the 
source  <leplli  is  h.  The  ray  starts  with  the  p  —  sinS,  /0  and  reaches  the  surface 


(32)  and  evaluating  (29)  yields  a  simple  step  response,  as  discussed  by  Chapman  (1976). 
Essentially,  (33)  has  a  simple  square-root  singularity  at  r  =  r(p),  and  rays  that  hit  the 
surface  near  the  receiver  dominate  the  behavior.  Since  p  varies  along  the  path,  we  must 
define  which  p  to  use  in  the  evaluation  of  (32).  The  proper  choice  is  the  starting  p  at 
the  source  as  outlined  in  the  previous  section.  Note  that  for  the  case  of  an  up-going 
direct  ray  the  two  methods  can  be  interpreted  in  a  similar  manner.  Only  one  ray  is 
involved  in  both,  and  applying  the  first- motion  approximation  of  (30)  yields  (32)  where 
the  extra  (2)  is  produced  by  the  double  valued  nature  of  expression  (33).  Thus,  the 
application  of  WKBJ  theory  to  the  locally  dipping  problem  appears  to  be  essentially  the 
same  as  for  the  uniform  layered  problem.  We  trace  the  ray  through  a  stack  of  layers 
down  to  the  turning  region,  turn  it  around  analytically,  and  follow  it  to  the  surface 
obeying  Snell’s  law.  The  special  treatment  at  the  turning  point  removes  the  nonlinear 
ray  parameter  effects  of  the  homogeneous  layered  parameterization  as  is  well  known. 
Essentially,  when  the  ray  reaches  critical  angle  at  a  given  interface,  we  back  up  to  the 
previous  reference  interface  and  compute  the  local  linear  velocity  gradient.  The  travel 
time  and  location  at  which  the  ray  recrosses  the  reference  interface  is  easily  calculated 
by  analytical  means.  Such  a  procedure  is  compatible  with  the  Langer  approximation 
which  is  the  basis  for  the  WKBJ  method,  see  Aki  and  Richards  (1980).  Earth  stretching 
becomes  slightly  more  complicated  than  in  the  uniform  layered  case.  In  this  situation  we 
let 


raj  =  (ry  -h  ry+i)/2 

where  ry  is  the  radius  of  the  Earth  at  the  j  th  layer,  etc.,  and  roy  becomes  the  radius  at 
the  layer’s  midpoint.  The  cross-section  is  then  constructed  in  terms  of  vertical  profiles 
of  velocity  and  thickness  vs.  ray  and  the  points  connected  by  linear  lines  as  displayed  in 


Fig.  5.  Next,  the  layer  thickness,  Thj ,  are  increased  by 

Ty  =  Thj{r,/raj) 


K* 


h*. 
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where  r,  is  the  radius  of  the  Earth.  Note  that  at  this  stage  the  horizontal  velocity  in 
each  layer  remains  constant.  The  velocity-depth  cross-sections  will  be  presented  in  this 
format.  The  velocity  and  density  approximations  are  determined  as  the  ray  encounters 
the  various  interfaces  with 

0=  VjirJraj) 

Pj  =  <^iir,/Taj) 

and  ray  is  determined  by  the  local  layer  thickness  and  position.  Thus,  the  velocity  is  no 
longer  constant  in  any  given  layer  but  depends  on  local  depth  correction.  The  C,  factor 
can  be  assumed  to  be  one  for  most  applications  of  gentle  dipping  structure,  as  discussed 
in  the  next  section  and  was  omitted  from  (32). 

The  approach  followed  here  is  similar  to  that  followed  by  Wiggins  (1976)  in  that 
the  main  justification  for  expression  (33)  is  that  it  yields  results  comparable  to  GRT.  A 
theoretical  justification  of  applying  WKBJ  to  laterally  varying  structure  is  given  by 
Chapman  and  Drummond  (1982).  See  Wesson  (1970). 

3.  Applications 

In  this  section  we  will  briefly  outline  possible  applications  of  these  approximate 
solutions  to  seismological  problems.  First,  the  direct  or  up-going  energy  problem  is  dis¬ 
cussed  when  motions  in  the  sloping  layers  of  a  sedimentary  basin  are  excited  by  a  line- 
source.  In  this  form  Finite  Difference  calculations  can  be  used  to  check  the  accuracy  of 
the  GRT  results.  Next,  the  point  source  excitation  of  Love  waves  is  considered  in  the 
presence  of  sloping  structure  followed  by  models  of  growing  Lithosphere.  Finally,  we 
construct  synthetics  for  laterally  varying  upper  mantle  models  and  confirm  the  useful¬ 
ness  of  WKBJ  at  long  periods. 

Local  seismograms 


One  of  many  complexities  associated  with  strong  motion  seismology  is  the  notice¬ 
ably  long  duration  of  high  frequency  P-waves  observed  in  sedimentary  basins.  These 
waves  are  generally  polarized  onto  the  vertical  component  due  to  the  strong  velocity  gra¬ 
dients  near  the  surface.  The  latter  portion  of  these  observed  motions  are  generally 
depleted  at  lower  frequency.  Thus,  one  might  conclude  that  there  are  propagational 
waveguides  that  preferentially  prolong  high  frequency  motions.  Non-planar  surface 
layering  appears  to  have  this  property.  This  calculation  will  be  done  with  SH-waves 
since  this  type  of  motion  is  studied  throughout  the  remainder  of  this  paper,  but  we 
would  expect  that  P-waves  would  behave  in  a  similar  manner. 

A  single  low  velocity  layer  which  grows  with  distance  between  the  source  and 
receiver  is  assumed  with  a  line  source  of  SH  motion  situated  at  a  depth  of  5.5  km.  The 
response  build-up  as  a  function  of  the  number  of  multiples  is  displayed  in  Fig.  6.  The 
square-root  singularity  indicated  in  expression  (13)  is  apparent  for  the  direct  arrival. 
Note  that  after  one  bounce  the  reflection  from  the  lower  interface  becomes  complex 
because  of  the  local  ray  parameter  effect  and  a  head  wave  and  post-critical  angle 
reflection  develops.  After  two  bounces,  the  time  separation  between  the  head  wave 
onset  and  reflection  times  becomes  less  and  the  reflected  spike  increases  in  strength. 
.\fter  many  bounces  the  ray  can  no  longer  reach  critical  angle  and  still  fit  into  the 
waveguide.  Thus,  (R  )"  becomes  small  since  the  reflection  coefficient  (R)  becomes  less 
than  one.  The  drop-off  in  amplitude  of  the  multiples  occurs  abruptly  at  this  time  on  the 
record 

The  corresponding  point  source  response  displayed  in  Fig.  7  can  be  obtained  from 
expression  (29).  Neglecting  the  C,  factor  produces  a  similar  response  with  a  slight 
reduction  in  later  arrivals,  roughly  13%  for  the  last  arrival.  Thus,  point  source  ampli¬ 
tudes  can  be  approximated  quite  well  by  scaling  line  source  results  by  the  square-root  of 
the  distance  factor  similar  to  the  flat  case.  Note  that  the  Cagniard-de  Hoop  technique 
proves  particularly  useful  in  tracing  these  rays  and  evaluating  their  individual 


Figure  6:  Line  source  response  as  a  function  of  rav  summation.  The  seismic  parame¬ 
ters  are  Pi  =  1.5  km/sec,  Pi  —  1.5  gm  /cm^,  and  /?2  =  3.3,  po  —  2.5  for  the 
upper  and  lower  layer  respectively.  The  top  plot  displays  the  ray  paths  at  ray 
parameters  appropriate  for  Snell’s  law  at  true  scale.  Ray  plots  on  the  left  have  a 
vertical  exaggeration  of  3  to  1. 
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coatributions.  However,  aa  mentioned  earlier,  this  series  of  rays  does  not  necessarily 
converge  to  the  exact  solution  and  some  demonstration  of  accuracy  is  required.  This 
was  attempted  earlier  by  Hong  and  Helmberger  (1977)  but  not  very  convincingly.  A 
much  more  rigorous  comparison  with  a  numerical  code  is  presently  being  conducted  by 
Vidale  et  al  (1984)  with  preliminary  results  of  the  comparison  of  the  two  techniques  for 
this  simple  model  displayed  in  Fig.  8.  The  top  trace  is  the  broadband  result  displayed  in 
Fig.  6.,  with  a  filtered  response  in  the  middle  for  comparison  with  finite-difference  results 
on  the  bottom.  The  highest  frequencies  have  been  removed  in  this  comparison  due  to 
computational  expense  but  the  existence  of  strong  high  frequency  multiples  is  striking. 
Since  the  finite-difference  calculation  can  be  performed  on  any  arbitrary  two-dimensional 
structure  we  have  extended  the  thin  layer  directly  above  the  source  to  the  left  as  a  flat 
thin  layer  avoiding  the  wedge  effect  which  is  obviously  omitted  in  the  ray  solution. 
Comparison  with  and  without  the  wedge  and  many  other  complexities  involving  double¬ 
couple  solutions  constructed  by  line-to-point  source  operators  are  discussed  in  Vidale  et 
al  (1984).  We  will  suppose  throughout  the  remainder  of  this  paper  that  the  generalized 
ray  modifications  discussed  in  the  previous  section  are  sufficiently  accurate  to  test  the 
WKBJ  synthetics. 


Another  interesting  application  of  the  above  technique  is  in  the  development  of 
Love  waves  and  the  effects  of  traveling  across  oceanic-to-continental  transitions.  This 
problem  was  encountered  in  a  recent  paper  by  Grand  and  Helmberger  (1981b)  when  the 
so-called  G-phase,  the  name  applied  to  the  impulsive  Love  waves  associated  with  oceanic 
paths,  interferes  with  mantle  arrivals.  Apparently,  this  situation  occurs  for  well 
developed  lithosphere  associated  with  older  plates  over-lying  slower  upper-mantle  velo¬ 
city  models.  The  beginning  portions  of  the  G-phase  as  recorded  slightly  inland  develop 
longer  periods  than  observed  at  island  stations.  Their  period  and  arrival  times  are 
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Figure  8:  Comparison  of  CRT  results  with  a  finite-difference  calculation.  The 
broadband  trace  has  been  filtered  to  remove  the  high  frequency  spikes. 
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compatible  with  the  model  presented  in  Fig.  9a.  A  dipping  model  with  arrival  times 
compatible  with  the  flat  model  is  presented  in  Fig.  9b  along  with  the  comparison  of  step 
responses  given  on  the  right.  Note  that  the  first  45  sec  of  motion  are  nearly  identical. 
The  higher  frequency  portions  of  the  Love  wave  become  less  pronounced  in  the  dipping 
case  but  the  general  appearance  is  similar  to  the  pure-oceanic  case,  see  Grand  and  Helm- 
berger  (1984b). 

It  appears  that  as  the  lithosphere  ages  it  gets  thicker,  for  example,  see  Forsyth 
(1975).  A  preliminary  model  of  predicted  Love  waves  for  this  situation  is  given  in  Fig. 
10,  also  included  are  synthetics  for  a  fast  and  slow  mantle.  The  long  period  nature  of 
the  synthetics  from  the  dipping  model  is  similar  to  the  slow  model  as  we  might  expect 
However,  there  is  considerable  roughness  at  the  start  of  the  Love  waves  caused  by  the 
mixed  paths  involving  both  the  crust  and  lid. 

Observationally,  we  see  upper-mantle  arrivals  starting  near  these  ranges.  Thus,  the 
diving  energy  must  be  added  to  these  synthetics  following  the  strategy  discussed  earlier. 
This  can  be  accomplished  by  summing  GR’s  or  by  applying  WKBJ. 

Upper -mantle  models 

In  this  section  we  investigate  effects  of  lateral  variation  in  upper  mantle  models,  as 
displayed  in  Fig.  11.  We  have  chosen  a  particularly  simple  case  with  no  low  velocity 
zone  to  simplify  the  comparison  of  GRT  with  WKBJ  synthetics.  A  further  simplification 
is  made  by  allowing  the  two  models  to  be  connected  in  a  linear  fashion  as  displayed  in 
the  middle  column. 

Following  the  WKBJ  approach  we  first  illuminate  the  model  by  tracking  a  set  of 
rays  from  the  source  towards  prospective  receivers.  These  rays  reach  the  .surface  at  r(p) 
in  time  T(p).  The  travel  time  at  a  particular  receiver,  r,  can  be  written  t(p)  ==-  p  r  -f 
T(p)  -  p  r  (p).  Note  that  p  changes  in  each  layer  but  they  are  all  functions  of  the  begin¬ 
ning  p.  Thus,  we  can  construct  the  t  versus  p  curves  as  displayed  in  Fig.  12  for  reversed 
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Figure  10:  Comparison  of  step  responses  for  two  simple  models  with  a  model  con¬ 
taining  a  growing  high  velocity  lid.  The  dotted  lines  indicate  the  first  arrival 
times  of  the  two  simple  models. 
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profiles.  The  largest  ray  parameter,  Pmwo  -26  which  corresponds  to  the  crustal  velo¬ 
city  of  3.9.  Next,  we  perform  the  numerical  derivative  (6p/St)  of  these  curves.  Note 
that  there  will  be  a  large  truncation  phase  at  the  near  stations  at  Pn»x-  This  can  be 
avoided  by  including  the  product  of  the  transmission  coefficients,  TC’s,  across  the  Moho, 
the  reference  interface  discussed  earlier,  since  TC(p„„)  is  small.  Thus,  the  product  of 
the  TC’s  with  (6p  jbt )  has  a  relatively  smooth  behavior.  The  head  wave  along  the  Moho 
is  added  in  by  including  the  reflected  generalized  ray.  By  performing  the  convolution 
indicated  in  expression  (29)  we  derive  step  responses  from  (t  vs.  p  )  curves  displayed  in 
Fig.  12.  These  results  are  shown  in  Fig.  13.  Short  period  synthetics  are  included  to 
emphasize  the  rapid  decay  of  amplitude  at  the  triplication  tips.  Eliminating  the  trunca¬ 
tion  phase  discussed  here  can  also  be  achieved  by  a  modification  of  the  Gaussian  beam 
technique  as  developed  by  Madariaga  and  Papadimitriou  (1984). 

The  synthetics  at  the  smallest  ranges  are  completely  controlled  by  the  shallow 
structure  and  the  local  model.  Thus,  the  first  arrival  from  the  Fast-to-Slow  synthetics 
has  a  shorter  travel  time  which  causes  the  triplication  from  the  400  km  discontinuity  to 
arrive  later  than  in  the  reverse  profile. 

A  more  detailed  plot  of  the  Slow-to-Fast  profile  is  displayed  in  Fig.  14  along  with 
the  GRT  responses  for  comparison.  The  synthetics  are  appropriate  for  the  WWSSN 
long  period  system.  A  typical  strike-slip  source  was  assumed  with  a  triangular  time  his¬ 
tory  of  (1,  1,  1  secs)  and  a  <*  =3,  see  Grand  and  Helmberger  (1984a). 

Note  that  there  is  a  distinct  change  in  the  latter  portion  of  the  WKBJ  step 
responses  between  17®  and  18® .  This  is  caused  by  omitting  the  head  waves  from 
along  the  top  of  the  model  for  distances  beyond  17*  .  However,  no  apparent  change  in 
the  synthetics  occurs  at  this  range  suggesting  that  the  long  period  drift  is  outside  the 
pass-band  of  the  operators  used  in  generating  these  .synthetics.  The  high-frequency 
spikes  so  apparent  in  the  GRT  step  responses  are  like- wise  removed  by  the  convolution 
operators. 
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Figure  14:  Step  responses  and  synthetics  generated  by  the  WKBJ  and  GRT  methods 
for  the  Slow-to-Fast  model  displayed  in  Fig.  11.  The  numbers  associated  with 
each  trace  indicate  the  maximum  amplitude. 
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The  400  km  discontinuity  is  treated  slightly  differently  in  the  two  methods  which 
leads  to  some  shifts  in  the  triplication  position.  In  GRT,  the  400  km  discontinuity  is 
generally  treated  as  a  sharp  jump  in  velocity  since  this  leads  to  the  best  results  when 
compared  against  reflectivity,  see  Burdick  and  Orcutt  (1978).  On  the  other  hand,  \VKBJ 
requires  a  smooth  transition,  3  km  transition  in  this  particular  case,  such  that  the  (p  vs. 
t  )  curve  is  smooth.  Thus,  the  sharp  spikes  occurring  in  the  GRT  step  responses  near 
14®  are  precritical  angle  reflections  from  the  400  km  discontinuity.  Similarly,  the  tripli¬ 
cation  seems  to  extend  to  greater  distances  in  the  GRT  results.  Note  that  the  most 
severe  mismatch  occurs  near  this  range.  At  larger  ranges  the  two  methods  agree  quite 
well,  especially  the  synthetic  waveforms.  In  fact,  the  synthetic  waveforms  agree  at  all 
distances  with  the  maximum  deviation  in  amplitude  of  about  25%.  And,  since  these 
synthetic  waveforms  are  used  to  interpret  observations  which  can  seldom  be  modeled  as 
well  as  the  agreement  between  these  two  methods,  we  can  consider  the  VVKBJ 
modifications  successful.  For  more  precision  involving  sharp  boundaries  we  suggest 
breaking  the  p  integration  into  a  combination  of  NVKBJ  for  the  smooth  portion  of  the 
model  and  a  generalized  ray  for  the  reflecting  interface,  for  example  .see  Given  (1984). 

Conclusions 

In  this  paper  we  presented  a  hybrid  procedure  of  generating  complete  seismograms 
in  laterally  varying  structure  by  applying  a  mixture  of  GRT  and  WKBJ  methods.  First, 
we  reviewed  the  modifications  of  GRT  required  for  dipping  structure  in  terms  of  local 
coordinates  and  ray  parameter  concepts  for  line  and  point  source  theory.  Solutions  cal¬ 
culated  by  this  approach  not  only  agree  with  geometric  results,  but  al.so  agree  with 
longer  period  motions  such  as  computed  with  finite-difference  methods.  Using  the 
correspondence  between  GRT  and  WKBJ  theory,  we  can  express  the  latter  in  relatively 
simple  form,  essentially  applying  a  square-root  of  distance  correction  to  line  source 
spreading.  Comparisons  between  GRT  and  WKBJ  synthetics  of  diving  energy  paths 


agree  reasonably  well.  Thus,  we  can  construct  nearly  complete  seismograms  with  a  com¬ 
bination  of  GRT  and  WTvBJ  with  the  former  used  to  handle  the  shallow  structure. 
Some  useful  demonstrations  of  the  methods  are  given  for  crustal  and  upper-mantle 


models. 
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ni.  A  Stable  Free- Surface  Boundary  Condition  for  2D  Elastic 
Finite-Difference  Wave  Simulation 

One  of  the  persistent  problems  in  finite-difference  solutions  of  the  elastic  wave 
equation  is  the  limited  stability  range  of  the  free-surface  boundary  condition,  and  its 
behavior  with  lateral  variations  in  velocity  and  density.  The  centered  difference  approxi¬ 
mation  (CDA)  presented  by  Alterman  and  Karal  (1968),  for  example,  remains  stable  only 
for  /?/o(  >  0.30,  where  and  a  are  the  shear  and  compressional  wave  velocities.  The 
one-sided  approximation  (OSA)  (Alterman  and  Rotenberg,  1969)  and  composed  approxi¬ 
mation  (Ilan  et  al.,  1975)  have  similar  restrictions.  The  revised  composed  approximation 
(CA)  of  Ilan  and  Loewenthal  (1976)  has  overcome  this  restriction,  but  apparently  cannot 
properly  handle  laterally  varying  media. 

In  the  first  part  of  this  paper  we  present  a  free-surface  boundary  approximation 
that  is  itself  stable  for  all  physical  0j(k  ratios  and  is  correct  for  laterally  varying  media. 
When  the  proposed  boundary  conditions  are  coupled  with  second-  and  fourth-order 
approximations  for  the  elastic  wave  equation,  the  overall  problem  is  stable  for  Pfot 
greater  than  0.01  and  ^/ot  greater  than  0.02,  respectively.  A  simple  numerical  test  of 
the  method  and  a  comparison  with  other  published  methods  is  given  in  the  second  sec¬ 
tion. 


FREE-SURFACE  BOUNDARY  CONDITIONS 

The  two-dimensional  free-surface  boundary  conditions  of  zero  tangential  and  nor¬ 
mal  stress  are 


du 

dz 
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(1) 
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du 

dx 
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(2) 


where  u  and  w  are  the  horizontal  and  vertical  displacements,  x  and  z  are  the  horizon¬ 
tal  and  vertical  spatial  coordinates,  and  qr  is  (  1  -  2o?l0^  ).  To  apply  finite-differences 
to  equations  (1)  and  (2),  an  extra  row  is  introduced  above  the  actual  free  surface.  The 


geometry  is  shown  in  Figure  1.  The  standard  finite  difference  approximations  to  the 
elastic  wave  equation  (Kelly  et  al.,  1976)  can  be  used  to  determine  the  solution  on  the 
interior  the  mesh,  up  to  and  including  row  1.  Previous  Tree-surface  boundary  condi¬ 
tions  cited  above  have  used  explicit  finite-difference  approximations  to  determine  the  row 
0.  This  approach  can  lead  to  stability  restrictions  in  terms  of  the  ratio  and  prob¬ 
lems  with  laterally  varying  media  due  to  the  difficulty  in  centering  the  normal  and 
tangential  derivatives  at  the  same  time  as  averaging  laterally  varying  media  parameters. 

The  method  proposed  here  uses  an  implicit  formulatbn  which  centers  both  the  nor¬ 
mal  and  tangential  derivatives  at  the  free  surface  which  is  halfway  between  row  0  and 
row  1.  The  scheme  is  similar  in  concept  to  the  Crank-Nicholson  method  for  the  diffusion 
equation  (e.g.,  Claierbout,  1976,  p.  185).  Applying  centered  second-order  differences  to 
equations  (1)  and  (2),  we  have 

Uo  -  j  B  Wo  =  u,  -H  j  B  w,  (3) 

Wo  -  r  B  uo  =  Wi  -1-  y  r  B  u,  (4) 

4  4 

where  the  subscripts  denote  the  solution  on  rows  0  and  1,  F  is  a  diagonal  matrix  which 
contains  the  values  of  across  the  surface,  and  B  is  a  bi-diagonal  matrix  with  sub-  and 
super-diagonals  equal  to  -1  and  1  respectively.  Note  that  the  centering  of  the  x  deriva¬ 
tives  is  achieved  by  averaging  estimates  on  rows  0  and  1. 

Equations  (3)  and  (4)  can  be  reduced  to  separate  systems  for  the  vector  unknowns 
Uq  and  Wq 

(I-.iBrB)u„-(I  +  iBrB)u,  +  iBw,  (5) 

(i-.irB»)wo=.(i  +  irB*)w,  +  irBu,  (6) 

The  matrices  on  the  left  side  of  equations  (5)  and  (6)  are  pentadiagonal  and  can  be 
solved  rapidly  by  an  algorithm  that  is  a  simple  extension  of  the  standard  tridiagonal 
solver  (e.g.  Claerbout,  1976,  p.  189).  The  vectors  on  the  right  side  can  be  computed 


from  displacements  on  row  1. 


To  show  the  stability  of  equations  (S)  and  (6)  as  extrapolation  operators,  we 
Fourier  transform  over  x  .  This  leads  to  the  system 


where  u  and  w  now  denote  the  Fourier  duals  of  u  and  w,  it  is  the  dual  of  x ,  and  h  is 
the  horizontal  mesh  spacing.  The  eigenvalues  of  the  matrix  in  equation  (7)  are  unity 
which  means  that  the  boundary  conditions  are  stable.  To  show  that  the  combined  prob¬ 
lem  of  the  boundary  conditions  and  the  interior  solution  is  stable  is  beyond  the  scope  of 
this  note.  In  numerical  tests,  the  combined  problem  remains  stable  for  /9/a  greater  than 
0.01  for  a  second-order  interior  method,  and  0/a  greater  than  0.02  for  a  fourth-order 
method.  We  suspect  that  the  stability  problem  lies  with  the  interior  solutions  rather 
than  the  boundary  conditions  for  0/a  of  0. 


EDGES  OF  THE  FREE-SURFACE 

The  boundary  conditions  derived  above  must  be  modified  for  the  extreme  two  edge 
elements  on  each  side  of  the  free  surface,  which  are  shown  as  open  circles  in  Figure  1. 
For  these  points,  we  apply  the  Bl  absorbing  boundary  conditions  of  Clayton  and 
Engquist  (1977).  For  the  component  u  on  the  left  side  the  boundary  conditions  in  equa¬ 
tion  (5)  are  modified  to  be 

(  1  -H  6)u0^  -h  (  1  -  S)u/  =  (  1  -  1  (8) 

(  1  -I-  5)tti*  -I-  (  1  -  =  (  1  -  +  (  1  +  %2‘"‘  (9) 

where  6=  a  At  /  h,  and  h  is  the  mesh  spacing  and  At  is  the  time  step.  Here 
(  Vq,  « |,  1(2)  three  elements  of  the  vector  Uq,  and  the  superscripts  t  and 

t-1  refer  to  the  present  and  previous  time  step.  Stability  of  equations  (8)  and  (9)  is 
independent  of  the  0/a  ratio.  These  equations  most  effectively  absorb  horizontally  trav¬ 
eling  P-waves.  Similar  equations  are  used  for  the  vertical  component  «; ,  except  that 
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6^  {  0  At  )  /  h  k  used  to  absorb  horizontally  traveling  S-waves.  Mirror  images  of 
these  conditions  are  used  at  the  right  edge  the  free  surface. 

NUMERICAL  EXAMPLES  OF  STABILITY  AND  ACCURACY 

The  accuracy  of  the  proposed  implicit  approximation  (lA)  is  shown  in  Figure  2. 
The  CA  and  COA,  which  are  not  shown,  are  about  as  accurate  as  the  lA  because  all 
three  are  accurate  to  second-order.  The  OSA,  which  is  also  not  shown,  is  mtly  accurate 
to  first  order  and  therefore  introduces  more  error  at  the  shorter  wavelengths.  Both  the 
direct  and  the  Rayleigh  wave  agree  with  the  analytic  calculation.  The  slight  difference 
in  sharpness  between  the  analytic  and  the  finite  difference  traces  is  due  to  grid  dispersion 
in  the  finite  difference  method.  The  small  perturbations  to  the  Rayleigh  wave  are  shear 
wave  noise  which  comes  from  the  source  region. 

The  stability  of  the  proposed  free  surface  boundary  conditions  is  illustrated  by  cal¬ 
culations  for  the  vertical  component  of  a  half-space  problem  where  01  a  k  0.2.  These 
results  are  shown  in  Figure  3.  The  CA  and  lA  are  well-behaved,  while  the  OSA  and 
CDA  are  not.  The  slight  difference  in  amplitude  between  the  CA  and  lA  results  from 
the  slight  difference  in  the  free  surface  position,  which  is  half  a  mesh  spacing  further 
from  the  source  with  the  CA  method  than  with  the  LA  method.  The  results  for  the  sta¬ 
bility  of  the  OSA,  CDA,  and  CA  are  in  i^reement  with  those  found  by  Ilan  and 
Loewenthal  (1976). 

A  test  with  lateral  heterogeneity  is  shown  in  Figure  4.  The  velocities  within  the 
rectangle  under  the  receiver  are  a  factor  of  three  less  than  the  velocities  in  the  halfspace. 
The  OSA  is  inaccurate,  particularly  on  the  horizontal  component.  The  CDA  is  unstable. 
The  CA  is  more  gently  unstable,  but  in  our  experience,  structure  with  more  lateral  vari¬ 
ation  causes  the  method  to  become  unstable  more  rapidly.  Part  of  the  disagreement 
between  the  CA  and  the  lA  in  the  early  part  of  the  record  arises  from  the  slightly 
different  location  of  the  free  surface  mentioned  above.  Only  the  lA  results  in  energy 
dying  away  with  time  in  the  signal.  We  suspect  the  result  is  accurate,  but  we  have  no 
method  with  which  it  may  be  conveniently  checked. 
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Figure  2.  The  analytic  result  for  a  halfspace  with  an  explosive  Dne 
source  Is  compared  to  a  fourth-order  finite  difference  solution  with  the 
implicit  free-surface  boundary  condition.  The  analytic  result  is  calculated 
by  the  method  of  de  Hoop  (1060).  In  the  halfspace,  a  is  3.5  km/sec  and  fi 
b  1.5  km/sec. 


s  ^  wv 


1  km 


Receiver 


free  surface 


□  SoTirce 


»  I  0.4  sec 

Figure  S.  The  vertical  component  for  the  one-sided  (OSA),  central- 
difference  (CDA),  composed  (CA),  and  implicit  (lA)  approximations  of  the 
free-surface  boundary  conditions  are  shown  for  the  small  fija  ratio  of  0.2. 
Traces  are  dotted  where  they  go  off-scale.  Peak  amplitudes  are  given  to 
the  right  of  each  trace.  In  the  haifspace,  a  is  3.5  km/sec  and  $  is  0.7 
km /sec.  An  explosive  line  source  is  used. 
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Figure  4.  The  horizontal  and  vertical  component  for  the  one-sided 
(OSA),  central-difference  (CDA),  composed  (CA),  and  implicit  (lA)  approx¬ 
imations  of  the  free-surface  boundary  conaitions  are  shown  for  laterally 
varying  structure.  In  the  hatched  region,  which  is  2.0  km  wide  and  0.1  km 
deep,  a  is  1.3  km/sec,  fi  is  0.6  km/sec,  and  p  is  1.0  g/cc.  In  the  rest  of  the 
usnace. 


jjjL 


-28- 


The  only  free  surface  which  is  stable  for  both  low  fifa  ratios  and  for  lateral  hetero¬ 
geneity  is  the  implicit  scheme  proposed  in  this  note. 
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